module mod_chem_thermo

   use iso_fortran_env, only: dp => real64
   use mod_database, only: TArray
   use stdlib_string_type
   use stdlib_optval
   implicit none
   private
   public :: type_chem_thermo
#include <chem.f90>
   real(dp), parameter :: ONE = 1.0_dp
   real(dp), parameter :: TWO = 2.0_dp
   real(dp), parameter :: THREE = 3.0_dp
   real(dp), parameter :: FOUR = 4.0_dp
   real(dp), parameter :: FIVE = 5.0_dp
   real(dp), parameter :: SIX = 6.0_dp
   real(dp), parameter :: SEVEN = 7.0_dp
   real(dp), parameter :: EIGHT = 8.0_dp
   real(dp), parameter :: NINE = 9.0_dp
   real(dp), parameter :: ONEOverTWO = ONE/TWO
   real(dp), parameter :: ONEOverTHREE = ONE/THREE
   real(dp), parameter :: ONEOverFOUR = ONE/FOUR
   real(dp), parameter :: ONEOverFIVE = ONE/FIVE

   type :: type_chem_thermo !! 基本单位都是 J, kmol, K
      character(len=:), allocatable :: name
      real(dp), allocatable :: temperature_ranges(:)
      real(dp), allocatable :: coefficients(:, :)
      real(dp) :: h0, delta_h0, T0
      real(dp) :: cp0_R, dcp_R

      real(dp) :: a(4)

   contains

      procedure :: initial => initial_thermo

      procedure :: get_cp_R => ThermoGetCpR
      procedure :: get_enthalpy_RT => ThermoGetEnthalpyRT
      procedure :: get_internal_energy_RT => ThermoGetInterEnergyRT
      procedure :: get_entropy_R => ThermoGetEntropyR
      procedure :: get_gibbs_RT => ThermoGetGibbsRT

   end type type_chem_thermo

contains

   subroutine initial_thermo(this, geometry)
      class(type_chem_thermo), intent(inout) :: this
      character(len=*), intent(in) :: geometry
      integer :: l

      real(dp) :: y0, dy0, y1, dy1

      if (geometry == "linear") then
         this%cp0_R = 3.5_dp
      elseif (geometry == "atom") then
         this%cp0_R = 2.5_dp
      else
         print *, "unknown geometry:", geometry
         error stop
      end if

      y0 = this%cp0_R
      dy0 = 0.0_dp
      y1 = this%get_cp_R(298.15_dp)
      dy1 = (this%get_cp_R(298.25_dp) - this%get_cp_R(298.15_dp))/ 0.1_dp

      this%a = GetFunCoeffs(y0, y1, dy0, dy1)

   end subroutine initial_thermo

   elemental function ThermoGetCpR(this, T) result(cp)
      class(type_chem_thermo), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: cp

      real(dp) :: T_, x

      integer :: n

      integer :: lc, l

      lc = 1
      T_ = min(T, this%temperature_ranges(size(this%temperature_ranges)))
      do l = 2, size(this%temperature_ranges) - 1
         if (T_ >= this%temperature_ranges(l)) lc = l
      end do

      select case (this%name)
      case ("NASA7")
         n = 7
         cp = sum(this%coefficients(1:5, lc)*Tarray(T_, n - 2))
      case ("NASA9")
         n = 9
         cp = sum(this%coefficients(1:7, lc)*(Tarray(T_, n - 2)/T_**2))
      end select

      if( T <= 250.0_dp) then
         cp = this%cp0_R
      elseif( T > 250.0_dp .and. T < 298.15_dp) then
         x = (T - 250.0_dp)/(298.15_dp - 250.0_dp)
         cp = this%a(4)*x**3 + this%a(3)*x**2 + this%a(2)*x + this%a(1)
      endif

   end function ThermoGetCpR

   elemental function ThermoGetInterEnergyRT(this, T) result(u)
      class(type_chem_thermo), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: u

      u = this%get_enthalpy_RT(T) - 1.0_dp

   end function ThermoGetInterEnergyRT

   elemental function ThermoGetEnthalpyRT(this, T, flg) result(h)
      class(type_chem_thermo), intent(in) :: this
      real(dp), intent(in) :: T
      logical, intent(in), optional :: flg
      real(dp) :: h
      real(dp) :: x

      real(dp) :: T_

      integer :: n

      integer :: lc, l

      lc = 1
      T_ = min(T, this%temperature_ranges(size(this%temperature_ranges)))

      do l = 2, size(this%temperature_ranges) - 1
         if (T_ >= this%temperature_ranges(l)) lc = l
      end do

      select case (this%name)
      case ("NASA7")
         n = 7
         h = sum(this%coefficients(1:6, lc)*[Tarray(T_, 5), ONE/T_]* &
                 [ONE, ONEOverTWO, ONEOverTHREE, ONEOverFOUR, ONEOverFIVE, ONE])
      case ("NASA9")
         n = 9
         h = sum(this%coefficients(1:8, lc)* &
                 [Tarray(T_, 7)/T_**2, ONE/T_]* &
                 [-ONE, log(T_), ONE, ONEOverTWO, ONEOverTHREE, ONEOverFOUR, ONEOverFIVE, ONE])
      end select

      if (T >= T_) h = h + this%get_cp_R(T_)*(T - T_)/T

   end function ThermoGetEnthalpyRT

   elemental function ThermoGetEntropyR(this, T) result(s)
      class(type_chem_thermo), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: s

      real(dp) :: T_
      integer :: n

      integer :: lc, l

      lc = 1
      T_ = min(T, this%temperature_ranges(size(this%temperature_ranges)))

      do l = 2, size(this%temperature_ranges) - 1
         if (T_ >= this%temperature_ranges(l)) lc = l
      end do

      select case (this%name)
      case ("NASA7")
         n = 7
         s = sum([this%coefficients(1:5, lc), this%coefficients(7, lc)]* &
                 [Tarray(T_, 5), ONE]* &
                 [log(T_), ONE, ONEOverTWO, ONEOverTHREE, ONEOverFOUR, ONE])
      case ("NASA9")
         n = 9
         s = sum([this%coefficients(1:7, lc), this%coefficients(9, lc)]* &
                 [TArray(T_, 7)/T_**2, ONE]* &
                 [-ONEOverTWO, -ONE, log(T_), ONE, ONEOverTWO, ONEOverTHREE, ONEOverFOUR, ONE])
         ! s = -this%coefficients(1, lc)*T_**(-2)*0.5_dp &
         !     -this%coefficients(2, lc)*T_**(-1) &
         !     +this%coefficients(3, lc)*log(T_) &
         !     +this%coefficients(4, lc)*T_ &
         !     +this%coefficients(5, lc)*T_**2*ONEOverTWO &
         !     +this%coefficients(6, lc)*T_**3*oneOverThree &
         !     +this%coefficients(7, lc)*T_**4*ONEOverFOUR &
         !     +this%coefficients(9, lc)
      end select

      if (T >= T_) s = s + this%get_cp_R(T_)*log(T/T_)

   end function ThermoGetEntropyR

   elemental function ThermoGetGibbsRT(this, T) result(g)
      class(type_chem_thermo), intent(in) :: this
      real(dp), intent(in) :: T
      real(dp) :: g

      g = this%get_enthalpy_RT(T, .false.) - this%get_entropy_R(T)

   end function ThermoGetGibbsRT


   function GetFunCoeffs(y0, y1, dy0, dy1) result(a)
      real(dp), intent(in) :: y0, y1, dy0, dy1
      real(dp) :: a(4)

      a(1) = y0
      a(2) = dy0
      a(3) = -3.0_dp*y0 - 2.0_dp*dy0 + 3.0_dp*y1 - dy1
      a(4) = 2.0_dp*y0 + dy0 - 2.0_dp*y1 + dy1

   end function

end module mod_chem_thermo
